Quant STM

library(here)
library(DBI)
library(RSQLite)
library(raster)
library(sf)
library(httr) # generic webservice package
library(tmap)

#NASIS tables
# pedon  to siteobs: siteobsiidref==siteobsiid

pedon <- data.table::fread("C:/R_Drive/Data_Files/LPKS_Data/Data/Soil_Pedon_Databases/NRCS/NASIS/NASIS_APRIL_2017/CSV_files/pedon.csv",
    sep = "|", header = T)

siteobs <- data.table::fread("C:/R_Drive/Data_Files/LPKS_Data/Data/Soil_Pedon_Databases/NRCS/NASIS/NASIS_APRIL_2017/CSV_files/siteobs.csv",
    sep = "|", header = T)


nasis_pedon <- pedon |> dplyr::select(peiid, siteobsiidref, earthcovkind1, earthcovkind2, pedlabsampnum) |> dplyr::left_join(siteobs |> dplyr::select(siteobsiid, obsdate), by=c('siteobsiidref'='siteobsiid'))


#join to LT_phy_chem

LT_phy_chem_nasis <- LT_phy_chem |> dplyr::left_join(nasis_pedon, by="peiid") |> dplyr::distinct()

tmap_mode("view")

tmap mode set to interactive viewing

  tm_shape(states) +
    tm_borders() +
    tm_lines(lwd = "strokelwd", legend.lwd.show = FALSE) +
    tm_shape(LT_data_points) +
    tm_dots(col="#1B9E77", size=0.3) +
    
    #tm_scale_bar(position = c(0.06, 0.05)) +
    tm_add_legend('fill', 
    col = c( "#1B9E77"),
    border.col = "grey40",
    size = 1,
    labels = c('KSSL'),
    title="Potential Sampling Areas") +
    tm_layout(main.title = "Loamy Tableland ESD", bg.color = "white", legend.outside = TRUE) +
    tm_view(set.view = c(-99.22, 39.13,  4))

LT_data_points_nlcd_data <- dplyr::bind_cols(LT_phy_chem_nasis, Lt_nlcd.df)

LT_data_points_nlcd <- LT_data_points_nlcd_data |> dplyr::filter(!is.na(texture_lab), !hzn_bot < hzn_top, !hzn_bot == hzn_top)

aqp::depths(LT_data_points_nlcd) <-  pedon_key ~ hzn_top  + hzn_bot 

converting profile IDs from integer to character

LT_data_points_nlcd.slab.r <- aqp::slab(LT_data_points_nlcd, fm = pedon_key ~ clay_total + silt_total + sand_total + bulk_density_third_bar + ph_h2o + organic_carbon_walkley_black + estimated_organic_carbon + aggregate_stability_05_2_mm, slab.structure=c(0,15), slab.fun=mean_na)

Note: aqp::slice() will be deprecated in aqp version 2.0
--> Please consider using the more efficient aqp::dice()

LT_data_points_nlcd.slab.r.1 <- LT_data_points_nlcd.slab.r |> dplyr::select(-c(contributing_fraction)) |> tidyr::pivot_wider(names_from=variable, values_from=value)
LT_data_points_nlcd.slab.r.1[ is.na(LT_data_points_nlcd.slab.r.1) ] <- NA
LT_data_points_nlcd.slab.r.1 <- LT_data_points_nlcd.slab.r.1 |> dplyr::rowwise() |> dplyr::mutate(soc = if_else(!is.na(organic_carbon_walkley_black), organic_carbon_walkley_black, estimated_organic_carbon))|> dplyr::mutate(soc_cl = soc_class(soc), ph_cl = ph_class(ph_h2o), db_cl = db_class(bulk_density_third_bar), txt_class = gettt(sand_total, silt_total, clay_total))


LT_data_points_complete <- LT_data_points_nlcd.slab.r.1 |> dplyr::filter(!is.na(soc_cl) & !is.na(db_cl) & !is.na(ph_cl) & !is.na(txt_class)) |> dplyr::select(pedon_key, top, bottom, soc_cl, ph_cl, db_cl, txt_class) |> dplyr::ungroup() |> dplyr::mutate(pedon_key = as.integer(pedon_key))

LT_data_points_nlcd_data_sub <- LT_data_points_nlcd_data |> dplyr::select(pedon_key, SSL_name, lat, lon, earthcovkind1,obsdate,lc) |> dplyr::distinct()

LT_data_points_complete <- LT_data_points_complete |> dplyr::left_join(LT_data_points_nlcd_data_sub, by="pedon_key")

#DSP4SH database

dsp_data_points_nlcd_data <- dplyr::bind_cols(pedon_lab, dsp_nlcd.df)

aqp::depths(dsp_data_points_nlcd_data) <-  pedon_ID ~ lay_depth_to_top  + lay_depth_to_bottom 

converting profile IDs from integer to character

dsp_data_points_nlcd.slab.r <- aqp::slab(dsp_data_points_nlcd_data, fm = pedon_ID ~ clay_tot_psa + silt_tot_psa + sand_tot_psa + Bulk_Density + ph_h2o + SOC_pct + KSSL_WSA, slab.structure=c(0,15), slab.fun=mean_na)

Note: aqp::slice() will be deprecated in aqp version 2.0
--> Please consider using the more efficient aqp::dice()
Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.Warning: Bad horizonation detected, first matching horizon selected. Use strict=TRUE to enforce QA/QC.

dsp_data_points_nlcd.slab.r.1 <- dsp_data_points_nlcd.slab.r |> dplyr::select(-c(contributing_fraction)) |> tidyr::pivot_wider(names_from=variable, values_from=value)
dsp_data_points_nlcd.slab.r.1[ is.na(dsp_data_points_nlcd.slab.r.1) ] <- NA
dsp_data_points_nlcd.slab.r.1 <- dsp_data_points_nlcd.slab.r.1 |> dplyr::rowwise() |> dplyr::mutate(soc_cl = soc_class(SOC_pct), ph_cl = ph_class(ph_h2o), db_cl = db_class(Bulk_Density), txt_class = gettt(sand_tot_psa, silt_tot_psa, clay_tot_psa))
dsp_data_points_nlcd.slab.r.1 <- dsp_data_points_nlcd.slab.r.1 |> dplyr::ungroup() |> dplyr::rename(pedon_key=pedon_ID)

dsp_data_points_complete <- dsp_data_points_nlcd.slab.r.1 |> dplyr::filter(!is.na(soc_cl) & !is.na(db_cl) & !is.na(ph_cl) & !is.na(txt_class)) |> dplyr::select(pedon_key, top, bottom, soc_cl, ph_cl, db_cl, txt_class) |> dplyr::ungroup() |> dplyr::mutate(pedon_key = as.integer(pedon_key))

dsp_data_points_nlcd <- dplyr::bind_cols(pedon_lab, dsp_nlcd.df)
dsp_data_points_nlcd_data_sub <- dsp_data_points_nlcd |> dplyr::select(pedon_key=pedon_ID, SSL_name=Soil, lat=pedon_y, lon=pedon_x, earthcovkind1=trt,lc) |> dplyr::distinct()
LT_data_points_nlcd_data_sub <- LT_data_points_nlcd_data |> dplyr::select(pedon_key, SSL_name, lat, lon, earthcovkind1,obsdate,lc) |> dplyr::distinct()
dsp_data_points_complete <- dsp_data_points_complete |> dplyr::left_join(dsp_data_points_nlcd_data_sub, by="pedon_key")

combine dataframes

stm_data <- bind_rows(dsp_data_points_complete |> droplevels(), LT_data_points_complete |> dplyr::select(-c(obsdate))|> droplevels())

Error in `bind_rows()`:
! Can't combine `..1$ph_cl` <ordered<10807>> and `..2$ph_cl` <ordered<f70bd>>.
Backtrace:
  1. dplyr::bind_rows(...)
  5. vctrs (local) `<fn>`()
  6. vctrs::vec_default_ptype2(...)
 10. vctrs::stop_incompatible_type(...)
 11. vctrs:::stop_incompatible(...)
 12. vctrs:::stop_vctrs(...)

# Select features
mod <- train(ode_cl_aic_mlr_fs, task = t)

Error in UseMethod("trainLearner") : 
  no applicable method for 'trainLearner' applied to an object of class "c('bnc', 'RLearnerClassif', 'RLearner', 'Learner')"

Pixel /Index Value Description

1 No Change 2 Change from or to Water 3 Change from or to any of the four Urban classes (open space; low, medium, and high intensity) 4 Change from Herbaceous Wetland to Woody Wetland, or vice versa 5 Change from or to Herbaceous Wetland 6 Change from Cultivated Crops to Hay / Pasture, or vice versa 7 Change from or to Cultivated Crops 8 Change from or to Hay / Pasture 9 Persistent Grassland and Shrubland change. This change index attempts to identify changes to persistent Grassland and Shrubland areas, and to separate them from transitional shrubland areas such as regenerating forests. 10 Change from or to Barren 11 Change from or to any of the three Forest classes (Evergreen, Deciduous, and Mixed) 12 Change from or to Woody Wetland

National Land Cover Database Class Legend and Description

Class Value Classification Description

Water

11 Open Water- areas of open water, generally with less than 25% cover of vegetation or soil.

12 Perennial Ice/Snow- areas characterized by a perennial cover of ice and/or snow, generally greater than 25% of total cover.

Developed

21 Developed, Open Space- areas with a mixture of some constructed materials, but mostly vegetation in the form of lawn grasses. Impervious surfaces account for less than 20% of total cover. These areas most commonly include large-lot single-family housing units, parks, golf courses, and vegetation planted in developed settings for recreation, erosion control, or aesthetic purposes.

22 Developed, Low Intensity- areas with a mixture of constructed materials and vegetation. Impervious surfaces account for 20% to 49% percent of total cover. These areas most commonly include single-family housing units.

23 Developed, Medium Intensity -areas with a mixture of constructed materials and vegetation. Impervious surfaces account for 50% to 79% of the total cover. These areas most commonly include single-family housing units.

24 Developed High Intensity-highly developed areas where people reside or work in high numbers. Examples include apartment complexes, row houses and commercial/industrial. Impervious surfaces account for 80% to 100% of the total cover.

Barren

31 Barren Land (Rock/Sand/Clay) - areas of bedrock, desert pavement, scarps, talus, slides, volcanic material, glacial debris, sand dunes, strip mines, gravel pits and other accumulations of earthen material. Generally, vegetation accounts for less than 15% of total cover.

Forest

41 Deciduous Forest- areas dominated by trees generally greater than 5 meters tall, and greater than 20% of total vegetation cover. More than 75% of the tree species shed foliage simultaneously in response to seasonal change.

42 Evergreen Forest- areas dominated by trees generally greater than 5 meters tall, and greater than 20% of total vegetation cover. More than 75% of the tree species maintain their leaves all year. Canopy is never without green foliage.

43 Mixed Forest- areas dominated by trees generally greater than 5 meters tall, and greater than 20% of total vegetation cover. Neither deciduous nor evergreen species are greater than 75% of total tree cover.

Shrubland

51 Dwarf Scrub- Alaska only areas dominated by shrubs less than 20 centimeters tall with shrub canopy typically greater than 20% of total vegetation. This type is often co-associated with grasses, sedges, herbs, and non-vascular vegetation.

52 Shrub/Scrub- areas dominated by shrubs; less than 5 meters tall with shrub canopy typically greater than 20% of total vegetation. This class includes true shrubs, young trees in an early successional stage or trees stunted from environmental conditions.

Herbaceous

71 Grassland/Herbaceous- areas dominated by gramanoid or herbaceous vegetation, generally greater than 80% of total vegetation. These areas are not subject to intensive management such as tilling, but can be utilized for grazing.

72 Sedge/Herbaceous- Alaska only areas dominated by sedges and forbs, generally greater than 80% of total vegetation. This type can occur with significant other grasses or other grass like plants, and includes sedge tundra, and sedge tussock tundra.

73 Lichens- Alaska only areas dominated by fruticose or foliose lichens generally greater than 80% of total vegetation.

74 Moss- Alaska only areas dominated by mosses, generally greater than 80% of total vegetation.

Planted/Cultivated

81 Pasture/Hay-areas of grasses, legumes, or grass-legume mixtures planted for livestock grazing or the production of seed or hay crops, typically on a perennial cycle. Pasture/hay vegetation accounts for greater than 20% of total vegetation.

82 Cultivated Crops -areas used for the production of annual crops, such as corn, soybeans, vegetables, tobacco, and cotton, and also perennial woody crops such as orchards and vineyards. Crop vegetation accounts for greater than 20% of total vegetation. This class also includes all land being actively tilled.

Wetlands

90 Woody Wetlands- areas where forest or shrubland vegetation accounts for greater than 20% of vegetative cover and the soil or substrate is periodically saturated with or covered with water.

Denomination ph range

Ultra acid < 3.5 Extremely acid 3.5–4.4 Very strongly acid 4.5–5.0 Strongly acid 5.1–5.5 Moderately acid 5.6–6.0 Slightly acid 6.1–6.5 Neutral 6.6–7.3 Slightly alkaline 7.4–7.8 Moderately alkaline 7.9–8.4 Strongly alkaline 8.5–9.0 Very strongly alkaline > 9.0

SOC classes (4 classes)

Db classes

ph_class <- function(ph){
if(is.na(ph)){
  class = NA
}else if(ph <= 3.5){
  class = 'Ultra acid'
}else if(ph > 3.5 & ph <= 4.5){
  class = 'Extremely acid'
}else if(ph > 4.5 & ph <= 5){
  class = 'Very strongly acid'
}else if(ph > 5 & ph < 5.6){
  class = 'Strongly acid'
}else if(ph >= 5.6 & ph <= 6){
  class = 'Moderately acid'
}else if(ph > 6 & ph <= 6.5){
  class = 'Slightly acid'
}else if(ph > 6.5 & ph <= 7.3){
  class = 'Neutral'
}else if(ph > 7.3 & ph <= 7.8){
  class = 'Slightly alkaline'
}else if(ph > 7.8 & ph <= 8.4){
  class = 'Moderately alkaline'
}else if(ph > 8.4 & ph <= 9){
  class = 'Strongly alkaline'
}else if(ph > 9 & ph <= 9){
  class = 'Very strongly alkaline'
}else {class = NA}
return(class)
}

soc_class <- function(soc){
if(is.na(soc)){
  class = NA
}else if(soc <= 1){
  class = 'low'
}else if(soc > 1 & soc <= 2){
  class = 'moderate'
}else if(soc > 2 & soc <= 3){
  class = 'high'
}else if(soc > 3){
  class = 'very high'
}else {class = NA}
return(class)
}

db_class <- function(Db){
if(is.na(Db)){
  class = NA
}else if(Db <= 1){
  class = 'very loose'
}else if(Db > 1 & Db <= 1.2){
  class = 'loose'
}else if(Db > 1.2 & Db <= 1.4){
  class = 'normal'
}else if(Db > 1.4 & Db <= 1.6){
  class = 'hard'
}else if(Db > 1.6){
  class = 'compact'
}
return(class)
}

gettt <- function(sand, silt, clay){
  if(is.na(sand) | is.na(silt) | is.na(clay)){
    x = NA
  } else if((silt + 1.5 * clay) < 15){
    x = "Sand"
  } else if((silt + 1.5 * clay) >= 15 & (silt + 2.0 * clay) < 30){
    x = "Loamy sand"
  } else if((clay >= 7) & (clay <= 20) & (sand > 52) & ((silt + 2.0 * clay) >= 30)){
    x = "Sandy loam"
  } else if((clay < 7) & (silt < 50) & ((silt + 2.0 * clay) >= 30)){
    x = "Sandy loam"
  } else if((clay >= 7) & (clay <= 27) & (silt >= 28) & (silt < 50) & (sand <= 52)){
    x = "Loam"
  } else if(((silt >= 50) & (clay >= 12) & (clay < 27)) | ((silt >= 50) & (silt < 80) & (clay < 12))){
    x = "Silt loam"
  } else if((silt >= 80) & (clay < 12)){
    x = "Silt"
  } else if((clay >= 20) & (clay < 35) & (silt < 28) & (sand > 45)){
    x = "Sandy clay loam"
  } else if((clay >= 27) & (clay < 40) & (sand > 20) & (sand <= 45)){
    x = "Clay loam"
  } else if((clay >= 27) & (clay < 40) & (sand <= 20)){
    x = "Silty clay loam"
  } else if((clay >= 35) & (sand >= 45)){
    x = "Sandy clay"
  } else if((clay >= 40) & (silt >= 40)){
    x = "Silty clay"
  } else if((clay >= 40) & (sand <= 45) & (silt < 40)){
    x = "Clay"
  }
  return(x)
}